********************************************************************************************************************************
******************* Replication Code for The Decline of Religion and Its Rise in Electoral Politics*****************************
************************************ John Huber and Ahmed Ezzledin Mohamed *****************************************************
********************************************************************************************************************************
cd "C:/Users/Ahmed/Dropbox/Replication_Files/Replication Files for CPS"
****** Import data 
set matsize 8000
set more off
use "wvs_relig.dta", clear

****************** Figure (1)
//Panel (a) 
histogram belief, percent xla(1/10, valuelabel noticks) barw(0.6) xtitle(God's importance)sch(s2mono)
graph export spirituality_histogram.png

//Panel (b)
histogram attend_1, percent xla(1/4, valuelabel noticks) barw(0.6) xtitle(Church attendance)sch(s2mono)
graph export hist_attendance.png

****************** Figure (2)
//The production of this plot involves several steps

//(1) Define and label the four religious groups
gen spiritual=(belief==10&belief~=.)
replace spiritual=. if belief==.

gen religrp=.
replace religrp=1 if (attender==0 & spiritual==0)
replace religrp=2 if (attender==1 & spiritual==0)
replace religrp=3 if (attender==0 & spiritual==1)
replace religrp=4 if (attender==1 & spiritual==1)

tab religrp, gen(relgrp)

label define religrp 1 " Non-Religious" 2 "Non-Spiritual Participants"  3 "Non-Participant Spirituals" 4 "Participant Spirituals", replace
label values religrp religrp

//(2) Estimate the average size of each group per country
collapse (mean)  relgrp1 relgrp2 relgrp3 relgrp4 [aw=S017], by(S025)cw

//(3) Generate the plot
label var relgrp1 "Non-religious"
label var relgrp2 "Non-believing participants"
label var relgrp3 "Non-participant believers"
label var relgrp4 "Believing participants"
graph hbox relgrp1 relgrp2 relgrp3 relgrp4, legend(off) sch(s2mono) showyvars
graph export groupsize.png

************** Figure (3)
use "wvs_relig.dta", clear 

qui reg  sdantigay sdbelief sdattendance  belief_attend  sdinc sdeduc employed  sdage sdage2 notlook  orthodox_i protestant_i catholic_i /// 
female    typecity1 typecity2   i.year i.ctry_survey_r [aw=S017] ,  vce(cluster ctry_survey_r)
est sto gay

qui reg   sdantiab  sdbelief sdattendance  belief_attend   sdinc sdeduc employed  sdage sdage2 notlook  orthodox_i protestant_i catholic_i /// 
female    typecity1 typecity2 i.year   i.ctry_survey_r [aw=S017] ,  vce(cluster ctry_survey_r)
est sto abort

qui reg   sdind_re  sdbelief sdattendance  belief_attend   sdinc sdeduc employed  sdage sdage2 notlook  orthodox_i protestant_i catholic_i /// 
female    typecity1 typecity2  i.year  i.ctry_survey_r [aw=S017] ,  vce(cluster ctry_survey_r)
est sto resp

qui reg   sdinequality  sdbelief sdattendance  belief_attend    sdinc sdeduc employed  sdage sdage2 notlook  orthodox_i protestant_i catholic_i /// 
female    typecity1 typecity2  i.year  i.ctry_survey_r [aw=S017] ,  vce(cluster ctry_survey_r)
est sto inequal

coefplot ( gay \ abort \  resp \ inequal ), ///
keep(sdbelief sdattendance belief_attend female sdeducation) grid(none) nolab aseq xline(0) byopts(rows(2)) coeflabels(sdattendance = "Attendance" sdbelief= "Belief"  belief_attend= "Belief*Attendance" female ="Female" sdeducation ="Education") eqlabels("Homosexuality" "Abortion"  ///
"Responsibility" "Inequality") mcolor("black")sch(s2mono) xtitle("Coefficient")

graph export attitudes.png

*************** Figure (4) 
reg  sdvote_loc c.sdbelief##c.sdattendance  sdincome_imp sdeduc employed  sdage sdage2 notlook  orthodox_i protestant_i catholic_i /// 
female  typecity1 typecity2 i.year i.ctry_survey_r  [aw=S017] ,  vce(cluster ctry_survey_r)
est sto vote_loc

reg  sdvote_nat c.sdbelief##c.sdattendance  sdincome_imp sdeduc employed  sdage sdage2 notlook  orthodox_i protestant_i catholic_i /// 
female  typecity1 typecity2 i.year i.ctry_survey_r   [aw=S017] ,  vce(cluster ctry_survey_r)
est sto vote_nat


coefplot ( vote_nat \ vote_loc ), ///
keep(sdbelief sdattendance c.sdbelief#c.sdattendance  sdeducation) grid(none) nolab aseq xline(0) byopts(rows(2)) coeflabels(sdattendance = "Attendance" sdbelief= "Belief"  c.sdbelief#c.sdattendance= "Belief*Attendance" sdeducation ="Education") eqlabels("Vote National"  ///
"Vote Local") mcolor("black") sch(s2mono) xtitle("Coefficient")

graph export vote_plot.png

************* Figure (5)
//Import country-level data with voting fractionalization measures
use votingfrac.dta, replace

//Plotting the relationships for the full data and its regional subsets
twoway(scatter sdPVP sdcorr_part_spirit2, ml(iso3c)ytitle(PVP)xtitle(Religious Cohesion)sch(s2mono)legend(off)title(All countries (98 surveys)))(lfit  sdPVP sdcorr_part_spirit2 )
graph export pvpall.png, replace

twoway(scatter sdPVP sdcorr_part_spirit2 if neoeurope==1, ml(iso3c)ytitle(PVP)xtitle(Religious Cohesion)sch(s2mono)legend(off)title(Neo-Europe (38 surveys)))(lfit  sdPVP sdcorr_part_spirit2 )
graph export pvp_neoeurope.png, replace

twoway(scatter sdPVP sdcorr_part_spirit2 if easteurope==1, ml(iso3c)ytitle(PVP)xtitle(Religious Cohesion)sch(s2mono)legend(off)title(East/Central Europe (21 surveys)))(lfit  sdPVP sdcorr_part_spirit2 )
graph export pvp_easteurope.png, replace

twoway(scatter sdPVP sdcorr_part_spirit2 if latinam==1, ml(iso3c)ytitle(PVP)xtitle(Religious Cohesion)sch(s2mono)legend(off)title(Latin America (32 surveys)))(lfit  sdPVP sdcorr_part_spirit2 )
graph export pvp_latinam.png, replace

************* Table (1)
reg sdPVP sdcorr_part_spirit2   sdreligfrac4 easteur latinam neoeurope africa, cluster(country) 
eststo m1

reg sdPVP sdcorr_part_spirit2   sdreligfrac4  sdlndm sdlngdp_l5   sdgini_ml5  sdover65 sdstatefund_ch sdherfrel00   catholic orthodox easteur latinam neoeurope africa i.year, cluster(country) 
eststo m2

reg sdPVP sdcorr_part_spirit2   sdreligfrac4     sdlndm    sdover65  sdherfrel00     easteur latinam neoeurope africa  i.year  , cluster(country) 
eststo m3

reg sdPVP sdcorr_part_spirit2   sdreligfrac4    sdlndm     sdover65  sdherfrel00          sdmean_bel   easteur latinam neoeurope africa  i.year , cluster(country) 
eststo m4

reg sdPVP sdcorr_part_spirit2   sdreligfrac4    sdlndm     sdover65  sdherfrel00          sdmean_attend  easteur latinam neoeurope africa  i.year , cluster(country)
eststo m5

esttab m1 m2 m3 m4 m5, tex cells(b(star fmt(3)) p(par fmt(3)))  style(tex) star(+ 0.10 * 0.05 ** 0.01 *** 0.001)  ///
		keep(  sdcorr_part_spirit2 sdlndm sdlngdp_l5   sdgini_ml5   sdover65 sdstatefund_ch sdherfrel00 catholic orthodox sdreligfrac4  easteurope latinam africa neoeurope sdmean_bel sdmean_attend) ///
		stats(r2 N, fmt(%9.3f %9.0g) labels("R$^2$" Obs)) ///
		varlabels( sdcorr_part_spirit2 "{\sc religious cohesion}" sdlndm "{\sc dm}" sdlngdp_l5 "{\sc ln(gdp)}" sdgini_ml5 "\sc{gini}" ///
	    sdover65 "{\sc over 65}" sdstatefund_ch "{\sc state funding}" sdherfrel00 "{\sc frac. religions sects}" sdreligfrac4 "{\sc frac. religious groups}" catholic "{\sc catholic}" orthodox "{\sc orthodox}"  ///
		easteurope "{\sc east/central europe}" latinam "{\sc latin america}" africa "{\sc africa}" neoeurope "{\sc neo-europe}" sdmean_bel "{\sc mean belief}" sdmean_attend "{\sc mean attendance}")


************* Table (2)
//Import the data
use "wvs_relig.dta", clear

//GDP Models
xtmixed   sdattendance c.sdbelief##c.sdlngdp_l5     sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year [w=S017]  || ctry_survey_r:, vce(cluster ctry_survey_r)
est sto pa1

xtmixed   sdattendance c.sdbelief##c.sdlngdp_l5     sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year if income_imp<=4  [w=S017] || ctry_survey_r:, vce(cluster ctry_survey_r)
est sto pa2

xtmixed   sdattendance c.sdbelief##c.sdlngdp_l5     sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year if income_imp>=7 [w=S017] || ctry_survey_r:, vce(cluster ctry_survey_r)
est sto pa3

// Gini Models
xtmixed   sdattendance c.sdbelief##c.sdgini_nl5 sdlngdp_l5   sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year [w=S017]  || ctry_survey_r:, vce(cluster ctry_survey_r)
est sto pa4

xtmixed   sdattendance c.sdbelief##c.sdgini_nl5 sdlngdp_l5   sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year if income_imp<=4 [w=S017]  || ctry_survey_r:, vce(cluster ctry_survey_r)
est sto pa5


xtmixed   sdattendance c.sdbelief##c.sdgini_nl5 sdlngdp_l5   sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year if income_imp>=7 [w=S017]  || ctry_survey_r:, vce(cluster ctry_survey_r)
est sto pa6

// Redistribution Models
xtmixed   sdattendance c.sdbelief##c.sdredis_l5 sdlngdp_l5   sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year [w=S017]  || ctry_survey_r:, vce(cluster ctry_survey_r)
est sto pa7

xtmixed   sdattendance c.sdbelief##c.sdredis_l5 sdlngdp_l5   sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year if income_imp<=4 [w=S017]  || ctry_survey_r:, vce(cluster ctry_survey_r)
est sto pa8

xtmixed   sdattendance c.sdbelief##c.sdredis_l5 sdlngdp_l5   sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year if income_imp>=7 [w=S017]  || ctry_survey_r:, vce(cluster ctry_survey_r)
est sto pa9

estout    pa1 pa2 pa3 pa4 pa5 pa6 pa7 pa8 pa9, ///
cells(b(star fmt(3)) se(par fmt(3) )) ///
keep(sdbelief c.sdbelief#c.sdlngdp_l5 sdlngdp_l5 c.sdbelief#c.sdgini_nl5 sdredis_l5 c.sdbelief#c.sdredis_l5 sdgini_nl5 sdincome_imp sdeducation female) ///
varlabels(  sdbelief "{\sc belief}" sdredis_l5 "{\sc redistribution}" c.sdbelief#c.sdlngdp_l5 "{\sc belief*(ln)gdp}"   sdgini_nl5 "{\sc gini}" sdlngdp_l5 "{\sc (ln)gdp}" employed "{\sc employed}" female "{\sc female}"  sdincome_imp "{\sc ncome (std)" sdeducation "education (std)" notlooking "Out Labor Force" sdage "Age(std)" sdage2 "Age-sq(std)" orthodox_i "Orthodox" protestant_i "Protestant" catholic_i "Catholic" sdbelief "Belief(std)" sdattendance   "Attendance(std)" typecity1 "Rural/town" typecity2 "Small City" c.sdbelief#c.sdredis_l5 "{\sc belief*redistribution}" c.sdbelief#c.sdgini_nl5 "{\sc belief*gini(net)") ///
style(tex) ///
stats(N, fmt(0))

**************** Figure (6)
xtmixed   attendance c.belief##c.redis_l5 sdlngdp_l5   sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year if income_imp<=4 [w=S017]  || ctry_survey_r:, vce(cluster ctry_survey_r)	
predict a_poor if e(sample)

xtmixed   attendance c.belief##c.redis_l5 sdlngdp_l5   sdincome_imp sdage sdage2 notlook employed female sdeduc orthodox_i protestant_i catholic_i  sdstatefund_ch typecity1 typecity2 sdorthodox sdcatholic i.year if income_imp>=6 [w=S017]  || ctry_survey_r:, vce(cluster ctry_survey_r)	
predict a_rich if e(sample)

gen lowredis=redis_l5<=.089
replace lowredis=. if redis_l5>.089&redis_l5<.38
label var lowredis "Eq 1 for both quarter, 0 for top quarter, missing otherwise"
cap drop pa_*			  
gen pa_rich_low=.	
gen pa_poor_low=.	
gen pa_rich_high=.	
gen pa_poor_high=.		  
forvalues i=1/10{
*Low redistribution countries
	*poor people in low redis countries
	sum a_poor if belief ==`i'&redis_l5<=.089
	replace pa_poor_low=r(mean) if belief ==`i'&redis_l5<=.089
	* rich people in low redis countries
	sum a_rich if belief ==`i'&redis_l5<=.089
	replace pa_rich_low=r(mean) if belief ==`i'&redis_l5<=.089
*High redistribution countries
	*poor people in low redis countries
	sum a_poor if belief ==`i'&redis_l5<=.089
	replace pa_poor_high=r(mean) if belief ==`i'&redis_l5>=.38
	* rich people in low redis countries
	sum a_rich if belief ==`i'&redis_l5<=.089
	replace pa_rich_high=r(mean) if belief ==`i'&redis_l5>=.38
}

drop if lowredis==.
	
collapse a_rich a_poor, by(belief lowredis)	
	
scatter a_poor belief if lowre==1,legend(order(1 "Low redistribution" 2 "High Redistribution")) xtitle("Belief")ytitle(Mean attendance)sch(s2mono)title(Poorer individuals)||scatter a_poor belief if lowr==0	
graph export a_poor.pdf

scatter a_rich belief if lowre==1,legend(order(1 "Low redistribution" 2 "High Redistribution")) xtitle("Belief")ytitle(Mean attendance)sch(s2mono)title(Richer individuals)||scatter a_rich belief if lowr==0	
graph export a_rich.pdf
	
